My background

Proposed schedule

R packages

## Required packages for this workshop
lib2install <- c("metaSEM", "semPlot", "readxl")

## Install them automatically if they are not available on your computer
for (i in lib2install) {
  if (!(i %in% rownames(installed.packages()))) install.packages(i)
}

Some resources

Getting help in using MASEM

What is structural equation modeling (SEM)?

“for their development of models, statistical procedures, and a computer algorithm (LISREL by Karl G. Joreskog; EQS by Peter M. Bentler) for Structural Equation Modeling (SEM) that changed the way in which inferences are made from observational data by permitting hypotheses derived from theory to be tested.”

Problems in primary research using SEM

Terms used in the literature

Brown and Stayman (1992)

Premack and Hunter (1988)

Norton et al. (2013)

Murayama and Elliot (2012)

Approaches to MASEM

Univariate approach (Viswesvaran & Ones, 1995)

One (sample) size does not fit all

Problems of treating a correlation matrix as a covariance matrix

GLS approach (Becker, 1992)

GLS approach (2)

TSSEM approach

Fixed-effects TSSEM

Fixed-effects TSSEM: Stage 1 analysis (1)

Fixed-effects TSSEM: Stage 1 analysis (2)

Fixed-effects TSSEM: Stage 2 analysis (1)

Fixed-effects TSSEM: Stage 2 analysis (2)

Fixed-effects TSSEM: Stage 2 analysis (3)

## Asymptotic sampling covariance matrix based on the harmonic mean (375):
##           x2_x1     x3_x1     x3_x2
## x2_x1 0.0010923 0.0004864 0.0004864
## x3_x1 0.0004864 0.0010923 0.0004864
## x3_x2 0.0004864 0.0004864 0.0010923
## Asymptotic sampling covariance matrix calculated in TSSEM approach:
##           x2_x1     x3_x1     x3_x2
## x2_x1 0.0020480 0.0005768 0.0004079
## x3_x1 0.0005768 0.0008192 0.0002580
## x3_x2 0.0004079 0.0002580 0.0004096

Fixed-effects TSSEM: Stage 2 analysis (4)

Subgroup analysis

Random-effects TSSEM

Random-effects TSSEM: Stage 1 analysis (1)

Random-effects TSSEM: Stage 1 analysis (2)

Random-effects TSSEM: Stage 2 analysis

Correction for unreliability

Illustrations with R

A higher-order CFA model for the Big Five model

The dataset

Reading the Excel file into R

## Load the library to read XLSX file
library(readxl)

## Read the study characteristics
study <- read_xlsx("Digman97.xlsx", sheet="Info")

head(study)
## # A tibble: 6 x 3
##   Study                               n Cluster     
##   <chr>                           <dbl> <chr>       
## 1 Digman 1 (1994)                   102 Children    
## 2 Digman 2 (1994)                   149 Children    
## 3 Digman 3 (1963c)                  334 Children    
## 4 Digman & Takemoto-Chock (1981b)   162 Children    
## 5 Graziano & Ward (1992)             91 Adolescents 
## 6 Yik & Bond (1993)                 656 Young adults
## Create an empty list to store the correlation matrices
Digman97.data <- list()
  
## Read 1 to 14 correlation matrices
for (i in 1:14) {
  
  ## Read each sheet and convert it into a matrix
  mat <- as.matrix(read_xlsx("Digman97.xlsx", sheet=paste0("Study ", i)))
  
  ## Add the row names
  rownames(mat) <- colnames(mat)
  
  ## Save it into a list
  Digman97.data[[i]] <- mat
}

## Add the names of the studies
names(Digman97.data) <- study$Study

## Show the first few studies
head(Digman97.data)
## $`Digman 1 (1994)`
##        A     C   ES     E    I
## A   1.00  0.62 0.41 -0.48 0.00
## C   0.62  1.00 0.59 -0.10 0.35
## ES  0.41  0.59 1.00  0.27 0.41
## E  -0.48 -0.10 0.27  1.00 0.37
## I   0.00  0.35 0.41  0.37 1.00
## 
## $`Digman 2 (1994)`
##        A    C   ES     E     I
## A   1.00 0.39 0.53 -0.30 -0.05
## C   0.39 1.00 0.59  0.07  0.44
## ES  0.53 0.59 1.00  0.09  0.22
## E  -0.30 0.07 0.09  1.00  0.45
## I  -0.05 0.44 0.22  0.45  1.00
## 
## $`Digman 3 (1963c)`
##       A     C   ES     E    I
## A  1.00  0.65 0.35  0.25 0.14
## C  0.65  1.00 0.37 -0.10 0.33
## ES 0.35  0.37 1.00  0.24 0.41
## E  0.25 -0.10 0.24  1.00 0.41
## I  0.14  0.33 0.41  0.41 1.00
## 
## $`Digman & Takemoto-Chock (1981b)`
##        A     C   ES     E     I
## A   1.00  0.65 0.70 -0.26 -0.03
## C   0.65  1.00 0.71 -0.16  0.24
## ES  0.70  0.71 1.00  0.01  0.11
## E  -0.26 -0.16 0.01  1.00  0.66
## I  -0.03  0.24 0.11  0.66  1.00
## 
## $`Graziano & Ward (1992)`
##       A    C   ES    E    I
## A  1.00 0.64 0.35 0.29 0.22
## C  0.64 1.00 0.27 0.16 0.22
## ES 0.35 0.27 1.00 0.32 0.36
## E  0.29 0.16 0.32 1.00 0.53
## I  0.22 0.22 0.36 0.53 1.00
## 
## $`Yik & Bond (1993)`
##       A    C   ES    E    I
## A  1.00 0.66 0.57 0.35 0.38
## C  0.66 1.00 0.45 0.20 0.31
## ES 0.57 0.45 1.00 0.49 0.31
## E  0.35 0.20 0.49 1.00 0.59
## I  0.38 0.31 0.31 0.59 1.00
## Extract the sample sizes
Digman97.n <- study$n
Digman97.n
##  [1]  102  149  334  162   91  656   70   70  277  227 1000  227   91 1040
## Extract the cluster
Digman97.cluster <- study$Cluster
Digman97.cluster
##  [1] "Children"      "Children"      "Children"      "Children"      "Adolescents"   "Young adults"  "Young adults"  "Young adults"  "Mature adults"
## [10] "Mature adults" "Mature adults" "Mature adults" "Mature adults" "Mature adults"

Proposed model

model <- "## Factor loadings
          ## Alpha is measured by A, C, and ES
          Alpha =~ A + C + ES
          ## Beta is measured by E and I
          Beta =~ E + I
          ## Factor correlation between Alpha and Beta
          Alpha ~~ Beta"

## Display the model
plot(model, color="yellow")

## Convert the lavaan syntax into a RAM model as the metaSEM only knows the RAM model
## It is important to ensure that the variables are arranged in A, C, ES, E, and I.
RAM <- lavaan2RAM(model, obs.variables=c("A","C","ES","E","I"),
                  A.notation="on", S.notation="with")
RAM
## $A
##       A   C   ES  E   I   Alpha         Beta       
## A     "0" "0" "0" "0" "0" "0*AonAlpha"  "0"        
## C     "0" "0" "0" "0" "0" "0*ConAlpha"  "0"        
## ES    "0" "0" "0" "0" "0" "0*ESonAlpha" "0"        
## E     "0" "0" "0" "0" "0" "0"           "0*EonBeta"
## I     "0" "0" "0" "0" "0" "0"           "0*IonBeta"
## Alpha "0" "0" "0" "0" "0" "0"           "0"        
## Beta  "0" "0" "0" "0" "0" "0"           "0"        
## 
## $S
##       A          C          ES           E          I          Alpha             Beta             
## A     "0*AwithA" "0"        "0"          "0"        "0"        "0"               "0"              
## C     "0"        "0*CwithC" "0"          "0"        "0"        "0"               "0"              
## ES    "0"        "0"        "0*ESwithES" "0"        "0"        "0"               "0"              
## E     "0"        "0"        "0"          "0*EwithE" "0"        "0"               "0"              
## I     "0"        "0"        "0"          "0"        "0*IwithI" "0"               "0"              
## Alpha "0"        "0"        "0"          "0"        "0"        "1"               "0*AlphawithBeta"
## Beta  "0"        "0"        "0"          "0"        "0"        "0*AlphawithBeta" "1"              
## 
## $F
##    A C ES E I Alpha Beta
## A  1 0  0 0 0     0    0
## C  0 1  0 0 0     0    0
## ES 0 0  1 0 0     0    0
## E  0 0  0 1 0     0    0
## I  0 0  0 0 1     0    0
## 
## $M
##   A C ES E I Alpha Beta
## 1 0 0  0 0 0     0    0

Fixed-effects TSSEM: Stage 1 Analysis

## method="FEM": fixed-effects TSSEM
fixed1 <- tssem1(Digman97.data, Digman97.n, method="FEM")

## summary of the findings
summary(fixed1)
## 
## Call:
## tssem1FEM(Cov = Cov, n = n, cor.analysis = cor.analysis, model.name = model.name, 
##     cluster = cluster, suppressWarnings = suppressWarnings, silent = silent, 
##     run = run)
## 
## Coefficients:
##        Estimate Std.Error z value  Pr(>|z|)    
## S[1,2] 0.363278  0.013368 27.1760 < 2.2e-16 ***
## S[1,3] 0.390375  0.012880 30.3076 < 2.2e-16 ***
## S[1,4] 0.103572  0.015048  6.8830 5.862e-12 ***
## S[1,5] 0.092286  0.015047  6.1331 8.621e-10 ***
## S[2,3] 0.416070  0.012519 33.2344 < 2.2e-16 ***
## S[2,4] 0.135148  0.014776  9.1464 < 2.2e-16 ***
## S[2,5] 0.141431  0.014866  9.5135 < 2.2e-16 ***
## S[3,4] 0.244459  0.014153 17.2724 < 2.2e-16 ***
## S[3,5] 0.138339  0.014834  9.3259 < 2.2e-16 ***
## S[4,5] 0.424566  0.012375 34.3070 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                      Value
## Sample size                      4496.0000
## Chi-square of target model       1505.4443
## DF of target model                130.0000
## p value of target model             0.0000
## Chi-square of independence model 4471.4242
## DF of independence model          140.0000
## RMSEA                               0.1815
## RMSEA lower 95% CI                  0.1736
## RMSEA upper 95% CI                  0.1901
## SRMR                                0.1621
## TLI                                 0.6580
## CFI                                 0.6824
## AIC                              1245.4443
## BIC                               412.0217
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## extract coefficients
coef(fixed1)
##             A         C        ES         E          I
## A  1.00000000 0.3632782 0.3903748 0.1035716 0.09228557
## C  0.36327824 1.0000000 0.4160695 0.1351477 0.14143058
## ES 0.39037483 0.4160695 1.0000000 0.2444593 0.13833895
## E  0.10357155 0.1351477 0.2444593 1.0000000 0.42456626
## I  0.09228557 0.1414306 0.1383390 0.4245663 1.00000000

Fixed-effects TSSEM: Stage 2 Analysis (1)

fixed2 <- tssem2(fixed1, Amatrix=RAM$A, Smatrix=RAM$S, Fmatrix=RAM$F)
summary(fixed2)
## 
## Call:
## wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj), 
##     n = sum(tssem1.obj$n), Amatrix = Amatrix, Smatrix = Smatrix, 
##     Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## AonAlpha      0.562800  0.015380 0.532655 0.592944  36.593 < 2.2e-16 ***
## ConAlpha      0.605217  0.015324 0.575183 0.635251  39.495 < 2.2e-16 ***
## EonBeta       0.781413  0.034244 0.714297 0.848529  22.819 < 2.2e-16 ***
## ESonAlpha     0.719201  0.015685 0.688458 0.749943  45.852 < 2.2e-16 ***
## IonBeta       0.551374  0.026031 0.500354 0.602394  21.181 < 2.2e-16 ***
## AlphawithBeta 0.362678  0.022384 0.318806 0.406550  16.203 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                4496.0000
## Chi-square of target model                   65.4526
## DF of target model                            4.0000
## p value of target model                       0.0000
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model           3112.7380
## DF of independence model                     10.0000
## RMSEA                                         0.0585
## RMSEA lower 95% CI                            0.0465
## RMSEA upper 95% CI                            0.0713
## SRMR                                          0.0284
## TLI                                           0.9505
## CFI                                           0.9802
## AIC                                          57.4526
## BIC                                          31.8088
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)

Fixed-effects TSSEM: Stage 2 Analysis (2)

plot(fixed2, color="green")

Fixed-effects TSSEM with cluster: Stage 1 Analysis

# Display the original study characteristic
table(Digman97.cluster)     
## Digman97.cluster
##   Adolescents      Children Mature adults  Young adults 
##             1             4             6             3
## Younger participants: "Children" and "Adolescents"
## Older participants: "Mature adults"
sample <- ifelse(Digman97.cluster %in% c("Children", "Adolescents"), 
                 yes="Younger participants", no="Older participants")

table(sample)
## sample
##   Older participants Younger participants 
##                    9                    5
## cluster: variable for the analysis with cluster
fixed1.cluster <- tssem1(Digman97.data, Digman97.n, method="FEM", cluster=sample)

summary(fixed1.cluster)
## $`Older participants`
## 
## Call:
## tssem1FEM(Cov = data.cluster[[i]], n = n.cluster[[i]], cor.analysis = cor.analysis, 
##     model.name = model.name, suppressWarnings = suppressWarnings)
## 
## Coefficients:
##        Estimate Std.Error z value  Pr(>|z|)    
## S[1,2] 0.297471  0.015436 19.2716 < 2.2e-16 ***
## S[1,3] 0.370248  0.014532 25.4774 < 2.2e-16 ***
## S[1,4] 0.137694  0.016403  8.3944 < 2.2e-16 ***
## S[1,5] 0.098061  0.016724  5.8637 4.528e-09 ***
## S[2,3] 0.393692  0.014146 27.8306 < 2.2e-16 ***
## S[2,4] 0.183045  0.016055 11.4009 < 2.2e-16 ***
## S[2,5] 0.092774  0.016643  5.5743 2.485e-08 ***
## S[3,4] 0.260753  0.015554 16.7645 < 2.2e-16 ***
## S[3,5] 0.096141  0.016573  5.8009 6.596e-09 ***
## S[4,5] 0.411756  0.013900 29.6225 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                      Value
## Sample size                      3658.0000
## Chi-square of target model        825.9826
## DF of target model                 80.0000
## p value of target model             0.0000
## Chi-square of independence model 3000.9731
## DF of independence model           90.0000
## RMSEA                               0.1515
## RMSEA lower 95% CI                  0.1424
## RMSEA upper 95% CI                  0.1611
## SRMR                                0.1459
## TLI                                 0.7117
## CFI                                 0.7437
## AIC                               665.9826
## BIC                               169.6088
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## 
## $`Younger participants`
## 
## Call:
## tssem1FEM(Cov = data.cluster[[i]], n = n.cluster[[i]], cor.analysis = cor.analysis, 
##     model.name = model.name, suppressWarnings = suppressWarnings)
## 
## Coefficients:
##         Estimate Std.Error z value  Pr(>|z|)    
## S[1,2]  0.604330  0.022125 27.3141 < 2.2e-16 ***
## S[1,3]  0.465536  0.027493 16.9327 < 2.2e-16 ***
## S[1,4] -0.031381  0.035940 -0.8732   0.38258    
## S[1,5]  0.061508  0.034547  1.7804   0.07500 .  
## S[2,3]  0.501432  0.026348 19.0311 < 2.2e-16 ***
## S[2,4] -0.060678  0.034557 -1.7559   0.07911 .  
## S[2,5]  0.320987  0.031064 10.3330 < 2.2e-16 ***
## S[3,4]  0.175437  0.033675  5.2097 1.891e-07 ***
## S[3,5]  0.305149  0.031586  9.6609 < 2.2e-16 ***
## S[4,5]  0.478640  0.026883 17.8045 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                      Value
## Sample size                       838.0000
## Chi-square of target model        346.2810
## DF of target model                 40.0000
## p value of target model             0.0000
## Chi-square of independence model 1470.4511
## DF of independence model           50.0000
## RMSEA                               0.2139
## RMSEA lower 95% CI                  0.1939
## RMSEA upper 95% CI                  0.2355
## SRMR                                0.1411
## TLI                                 0.7305
## CFI                                 0.7844
## AIC                               266.2810
## BIC                                77.0402
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)

Fixed-effects TSSEM with cluster: Stage 2 Analysis (1)

fixed2.cluster <- tssem2(fixed1.cluster, Amatrix=RAM$A, Smatrix=RAM$S, Fmatrix=RAM$F)

summary(fixed2.cluster)
## $`Older participants`
## 
## Call:
## wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj), 
##     n = sum(tssem1.obj$n), Amatrix = Amatrix, Smatrix = Smatrix, 
##     Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## AonAlpha      0.512656  0.018206 0.476973 0.548340  28.158 < 2.2e-16 ***
## ConAlpha      0.549967  0.017945 0.514795 0.585140  30.647 < 2.2e-16 ***
## EonBeta       0.967136  0.058751 0.851986 1.082287  16.462 < 2.2e-16 ***
## ESonAlpha     0.732174  0.018929 0.695074 0.769273  38.680 < 2.2e-16 ***
## IonBeta       0.430649  0.029634 0.372568 0.488730  14.532 < 2.2e-16 ***
## AlphawithBeta 0.349236  0.028118 0.294125 0.404346  12.420 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                3658.0000
## Chi-square of target model                   21.9954
## DF of target model                            4.0000
## p value of target model                       0.0002
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model           2273.3236
## DF of independence model                     10.0000
## RMSEA                                         0.0351
## RMSEA lower 95% CI                            0.0217
## RMSEA upper 95% CI                            0.0500
## SRMR                                          0.0160
## TLI                                           0.9801
## CFI                                           0.9920
## AIC                                          13.9954
## BIC                                         -10.8233
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## 
## $`Younger participants`
## 
## Call:
## wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj), 
##     n = sum(tssem1.obj$n), Amatrix = Amatrix, Smatrix = Smatrix, 
##     Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##                Estimate Std.Error    lbound    ubound z value Pr(>|z|)    
## AonAlpha       0.747647  0.023880  0.700842  0.794451 31.3081   <2e-16 ***
## ConAlpha       0.911705  0.019864  0.872772  0.950638 45.8968   <2e-16 ***
## EonBeta        0.152561  0.159140 -0.159347  0.464469  0.9587   0.3377    
## ESonAlpha      0.677434  0.025864  0.626742  0.728126 26.1924   <2e-16 ***
## IonBeta        3.283875  3.363587 -3.308633  9.876384  0.9763   0.3289    
## AlphawithBeta  0.117255  0.125430 -0.128584  0.363094  0.9348   0.3499    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                 838.0000
## Chi-square of target model                  145.6167
## DF of target model                            4.0000
## p value of target model                       0.0000
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model           2480.2352
## DF of independence model                     10.0000
## RMSEA                                         0.2057
## RMSEA lower 95% CI                            0.1778
## RMSEA upper 95% CI                            0.2350
## SRMR                                          0.1051
## TLI                                           0.8567
## CFI                                           0.9427
## AIC                                         137.6167
## BIC                                         118.6926
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)

Fixed-effects TSSEM with cluster: Stage 2 Analysis (2)

## Setup two plots
layout(t(1:2))
invisible(lapply(fixed2.cluster, plot))

## Setup two plots
layout(t(1:2))

## Plot the first group
plot(fixed2.cluster[[1]])
title("Younger participants")

## Plot the second group
plot(fixed2.cluster[[2]])
title("Older participants")

Random-effects TSSEM: Stage 1 Analysis

## method="REM": Random-effects model
random1 <- tssem1(Digman97.data, Digman97.n, method="REM", RE.type="Diag")

summary(random1)
## 
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues, 
##     "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound, 
##     I2 = I2, model.name = model.name, suppressWarnings = TRUE, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##                Estimate   Std.Error      lbound      ubound z value  Pr(>|z|)    
## Intercept1   0.38971902  0.05429255  0.28330757  0.49613047  7.1781 7.068e-13 ***
## Intercept2   0.43265876  0.04145135  0.35141562  0.51390190 10.4377 < 2.2e-16 ***
## Intercept3   0.04945623  0.06071078 -0.06953471  0.16844718  0.8146   0.41529    
## Intercept4   0.09603703  0.04478710  0.00825593  0.18381813  2.1443   0.03201 *  
## Intercept5   0.42724236  0.03911620  0.35057603  0.50390870 10.9224 < 2.2e-16 ***
## Intercept6   0.11929312  0.04106203  0.03881303  0.19977322  2.9052   0.00367 ** 
## Intercept7   0.19292424  0.04757961  0.09966991  0.28617856  4.0548 5.018e-05 ***
## Intercept8   0.22690158  0.03240893  0.16338125  0.29042191  7.0012 2.538e-12 ***
## Intercept9   0.18105561  0.04258855  0.09758360  0.26452763  4.2513 2.126e-05 ***
## Intercept10  0.43614966  0.03205960  0.37331400  0.49898532 13.6043 < 2.2e-16 ***
## Tau2_1_1     0.03648370  0.01513054  0.00682839  0.06613902  2.4113   0.01590 *  
## Tau2_2_2     0.01963096  0.00859788  0.00277942  0.03648250  2.2832   0.02242 *  
## Tau2_3_3     0.04571436  0.01952284  0.00745029  0.08397842  2.3416   0.01920 *  
## Tau2_4_4     0.02236119  0.00995082  0.00285794  0.04186445  2.2472   0.02463 *  
## Tau2_5_5     0.01729071  0.00796404  0.00168149  0.03289994  2.1711   0.02992 *  
## Tau2_6_6     0.01815481  0.00895896  0.00059557  0.03571405  2.0264   0.04272 *  
## Tau2_7_7     0.02604881  0.01130265  0.00389601  0.04820160  2.3047   0.02119 *  
## Tau2_8_8     0.00988761  0.00513713 -0.00018098  0.01995619  1.9247   0.05426 .  
## Tau2_9_9     0.01988243  0.00895052  0.00233973  0.03742513  2.2214   0.02633 *  
## Tau2_10_10   0.01049221  0.00501578  0.00066147  0.02032295  2.0918   0.03645 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 1220.332
## Degrees of freedom of the Q statistic: 130
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                               Estimate
## Intercept1: I2 (Q statistic)    0.9326
## Intercept2: I2 (Q statistic)    0.8872
## Intercept3: I2 (Q statistic)    0.9325
## Intercept4: I2 (Q statistic)    0.8703
## Intercept5: I2 (Q statistic)    0.8797
## Intercept6: I2 (Q statistic)    0.8480
## Intercept7: I2 (Q statistic)    0.8887
## Intercept8: I2 (Q statistic)    0.7669
## Intercept9: I2 (Q statistic)    0.8590
## Intercept10: I2 (Q statistic)   0.8193
## 
## Number of studies (or clusters): 14
## Number of observed statistics: 140
## Number of estimated parameters: 20
## Degrees of freedom: 120
## -2 log likelihood: -112.4196 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Extract the fixed-effects estimates
(est_fixed <- coef(random1, select="fixed"))
##  Intercept1  Intercept2  Intercept3  Intercept4  Intercept5  Intercept6  Intercept7  Intercept8  Intercept9 Intercept10 
##  0.38971902  0.43265876  0.04945623  0.09603703  0.42724236  0.11929312  0.19292424  0.22690158  0.18105561  0.43614966
## Convert the estimated vector to a symmetrical matrix
## where the diagonals are fixed at 1 (for a correlation matrix)
vec2symMat(est_fixed, diag=FALSE)
##            [,1]      [,2]      [,3]       [,4]       [,5]
## [1,] 1.00000000 0.3897190 0.4326588 0.04945623 0.09603703
## [2,] 0.38971902 1.0000000 0.4272424 0.11929312 0.19292424
## [3,] 0.43265876 0.4272424 1.0000000 0.22690158 0.18105561
## [4,] 0.04945623 0.1192931 0.2269016 1.00000000 0.43614966
## [5,] 0.09603703 0.1929242 0.1810556 0.43614966 1.00000000

Random-effects TSSEM: Stage 2 Analysis

random2 <- tssem2(random1, Amatrix=RAM$A, Smatrix=RAM$S, Fmatrix=RAM$F)
summary(random2)
## 
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, Amatrix = Amatrix, 
##     Smatrix = Smatrix, Fmatrix = Fmatrix, diag.constraints = diag.constraints, 
##     cor.analysis = cor.analysis, intervals.type = intervals.type, 
##     mx.algebras = mx.algebras, model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## AonAlpha      0.569435  0.052425 0.466684 0.672187 10.8619 < 2.2e-16 ***
## ConAlpha      0.590630  0.052649 0.487439 0.693821 11.2182 < 2.2e-16 ***
## EonBeta       0.679955  0.075723 0.531541 0.828370  8.9795 < 2.2e-16 ***
## ESonAlpha     0.760455  0.061963 0.639009 0.881901 12.2726 < 2.2e-16 ***
## IonBeta       0.641842  0.072459 0.499825 0.783859  8.8580 < 2.2e-16 ***
## AlphawithBeta 0.377691  0.047402 0.284785 0.470596  7.9679 1.554e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                4496.0000
## Chi-square of target model                    7.8204
## DF of target model                            4.0000
## p value of target model                       0.0984
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model            501.6767
## DF of independence model                     10.0000
## RMSEA                                         0.0146
## RMSEA lower 95% CI                            0.0000
## RMSEA upper 95% CI                            0.0297
## SRMR                                          0.0436
## TLI                                           0.9806
## CFI                                           0.9922
## AIC                                          -0.1796
## BIC                                         -25.8234
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)

Plotting the results

## Plot the parameter estimates
plot(random2, color="green")

A regression model on SAT (Math) (1)

## Read the first sheet in the Excel file
my.xlsx <- read_excel("Becker94.xlsx")
my.xlsx
## # A tibble: 10 x 6
##    Study                    r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal     n gender 
##    <chr>                                 <dbl>                 <dbl>              <dbl> <dbl> <chr>  
##  1 Becker (1978) Females                 0.47                  -0.21              -0.15    74 Females
##  2 Becker (1978) Males                   0.28                   0.19               0.18   153 Males  
##  3 Berry (1957) Females                  0.48                   0.41               0.26    48 Females
##  4 Berry (1957) Males                    0.37                   0.4                0.27    55 Males  
##  5 Rosenberg (1981) Females              0.42                   0.48               0.23    51 Females
##  6 Rosenberg (1981) Males                0.41                   0.74               0.44    18 Males  
##  7 Weiner A (1984) Females               0.26                   0.72               0.36    27 Females
##  8 Weiner A (1984) Males                 0.32                   0.52               0.1     43 Males  
##  9 Weiner B (1984) Females               0.580                  0.64               0.4     35 Females
## 10 Weiner B (1984) Males                 0.34                   0.28              -0.03    34 Males
## Split the data by rows
my.R2 <- split(my.xlsx[, 2:4], seq(nrow(my.xlsx)))
head(my.R2)
## $`1`
## # A tibble: 1 x 3
##   r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
##                <dbl>                 <dbl>              <dbl>
## 1               0.47                 -0.21              -0.15
## 
## $`2`
## # A tibble: 1 x 3
##   r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
##                <dbl>                 <dbl>              <dbl>
## 1               0.28                  0.19               0.18
## 
## $`3`
## # A tibble: 1 x 3
##   r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
##                <dbl>                 <dbl>              <dbl>
## 1               0.48                  0.41               0.26
## 
## $`4`
## # A tibble: 1 x 3
##   r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
##                <dbl>                 <dbl>              <dbl>
## 1               0.37                   0.4               0.27
## 
## $`5`
## # A tibble: 1 x 3
##   r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
##                <dbl>                 <dbl>              <dbl>
## 1               0.42                  0.48               0.23
## 
## $`6`
## # A tibble: 1 x 3
##   r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
##                <dbl>                 <dbl>              <dbl>
## 1               0.41                  0.74               0.44
## Convert the row correlations into correlation matrices
my.R2 <- lapply(my.R2, function(x) vec2symMat(unlist(x), diag=FALSE))
head(my.R2)
## $`1`
##       [,1]  [,2]  [,3]
## [1,]  1.00  0.47 -0.21
## [2,]  0.47  1.00 -0.15
## [3,] -0.21 -0.15  1.00
## 
## $`2`
##      [,1] [,2] [,3]
## [1,] 1.00 0.28 0.19
## [2,] 0.28 1.00 0.18
## [3,] 0.19 0.18 1.00
## 
## $`3`
##      [,1] [,2] [,3]
## [1,] 1.00 0.48 0.41
## [2,] 0.48 1.00 0.26
## [3,] 0.41 0.26 1.00
## 
## $`4`
##      [,1] [,2] [,3]
## [1,] 1.00 0.37 0.40
## [2,] 0.37 1.00 0.27
## [3,] 0.40 0.27 1.00
## 
## $`5`
##      [,1] [,2] [,3]
## [1,] 1.00 0.42 0.48
## [2,] 0.42 1.00 0.23
## [3,] 0.48 0.23 1.00
## 
## $`6`
##      [,1] [,2] [,3]
## [1,] 1.00 0.41 0.74
## [2,] 0.41 1.00 0.44
## [3,] 0.74 0.44 1.00
## Add the labels
var.names <- c("SAT_Math", "Spatial", "SAT_Verbal")
my.R2 <- lapply( my.R2, function(x) { dimnames(x) <- list(var.names, var.names); x} )

## Add the study name
names(my.R2) <- my.xlsx$Study
head(my.R2)
## $`Becker (1978) Females`
##            SAT_Math Spatial SAT_Verbal
## SAT_Math       1.00    0.47      -0.21
## Spatial        0.47    1.00      -0.15
## SAT_Verbal    -0.21   -0.15       1.00
## 
## $`Becker (1978) Males`
##            SAT_Math Spatial SAT_Verbal
## SAT_Math       1.00    0.28       0.19
## Spatial        0.28    1.00       0.18
## SAT_Verbal     0.19    0.18       1.00
## 
## $`Berry (1957) Females`
##            SAT_Math Spatial SAT_Verbal
## SAT_Math       1.00    0.48       0.41
## Spatial        0.48    1.00       0.26
## SAT_Verbal     0.41    0.26       1.00
## 
## $`Berry (1957) Males`
##            SAT_Math Spatial SAT_Verbal
## SAT_Math       1.00    0.37       0.40
## Spatial        0.37    1.00       0.27
## SAT_Verbal     0.40    0.27       1.00
## 
## $`Rosenberg (1981) Females`
##            SAT_Math Spatial SAT_Verbal
## SAT_Math       1.00    0.42       0.48
## Spatial        0.42    1.00       0.23
## SAT_Verbal     0.48    0.23       1.00
## 
## $`Rosenberg (1981) Males`
##            SAT_Math Spatial SAT_Verbal
## SAT_Math       1.00    0.41       0.74
## Spatial        0.41    1.00       0.44
## SAT_Verbal     0.74    0.44       1.00

A regression model on SAT (Math) (2)

## Regression model
model <- "SAT_Math ~ Spatial + SAT_Verbal
          ## Correlation between Spatial and Verbal
          Spatial ~~ SAT_Verbal
          ## Fix the variances of the independent variables at 1.0
          Spatial ~~ 1*Spatial
          SAT_Verbal ~~ 1*SAT_Verbal"

## Plot the model
plot(model, color="yellow")

## Convert the equation into RAM 
RAM <- lavaan2RAM(model, obs.variables=c("SAT_Math", "Spatial", "SAT_Verbal"))
RAM
## $A
##            SAT_Math Spatial               SAT_Verbal              
## SAT_Math   "0"      "0*SAT_MathONSpatial" "0*SAT_MathONSAT_Verbal"
## Spatial    "0"      "0"                   "0"                     
## SAT_Verbal "0"      "0"                   "0"                     
## 
## $S
##            SAT_Math                 Spatial                   SAT_Verbal               
## SAT_Math   "0*SAT_MathWITHSAT_Math" "0"                       "0"                      
## Spatial    "0"                      "1"                       "0*SpatialWITHSAT_Verbal"
## SAT_Verbal "0"                      "0*SpatialWITHSAT_Verbal" "1"                      
## 
## $F
##            SAT_Math Spatial SAT_Verbal
## SAT_Math          1       0          0
## Spatial           0       1          0
## SAT_Verbal        0       0          1
## 
## $M
##   SAT_Math Spatial SAT_Verbal
## 1        0       0          0

Fixed-effects TSSEM: Stage 1 Analysis

## method="FEM": fixed-effects TSSEM
fixed1 <- tssem1(Becker94$data, Becker94$n, method="FEM")

## summary of the findings
summary(fixed1)
## 
## Call:
## tssem1FEM(Cov = Cov, n = n, cor.analysis = cor.analysis, model.name = model.name, 
##     cluster = cluster, suppressWarnings = suppressWarnings, silent = silent, 
##     run = run)
## 
## Coefficients:
##        Estimate Std.Error z value  Pr(>|z|)    
## S[1,2] 0.379961  0.037123 10.2351 < 2.2e-16 ***
## S[1,3] 0.334562  0.039947  8.3751 < 2.2e-16 ***
## S[2,3] 0.176461  0.042334  4.1683 3.069e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                      Value
## Sample size                       538.0000
## Chi-square of target model         63.6553
## DF of target model                 27.0000
## p value of target model             0.0001
## Chi-square of independence model  207.7894
## DF of independence model           30.0000
## RMSEA                               0.1590
## RMSEA lower 95% CI                  0.1096
## RMSEA upper 95% CI                  0.2117
## SRMR                                0.1586
## TLI                                 0.7709
## CFI                                 0.7938
## AIC                                 9.6553
## BIC                              -106.1169
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)

Random-effects TSSEM: Stage 1 Analysis

## method="REM": Random-effects model
random1 <- tssem1(Becker94$data, Becker94$n, method="REM", RE.type="Diag")
summary(random1)
## 
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues, 
##     "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound, 
##     I2 = I2, model.name = model.name, suppressWarnings = TRUE, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##              Estimate  Std.Error     lbound     ubound z value  Pr(>|z|)    
## Intercept1  0.3777491  0.0395030  0.3003246  0.4551736  9.5625 < 2.2e-16 ***
## Intercept2  0.3807843  0.0784956  0.2269357  0.5346329  4.8510 1.228e-06 ***
## Intercept3  0.1704927  0.0513545  0.0698398  0.2711457  3.3199 0.0009004 ***
## Tau2_1_1    0.0005038  0.0042009 -0.0077298  0.0087374  0.1199 0.9045407    
## Tau2_2_2    0.0416264  0.0257388 -0.0088206  0.0920734  1.6173 0.1058209    
## Tau2_3_3    0.0067540  0.0102792 -0.0133928  0.0269008  0.6571 0.5111466    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 61.02636
## Degrees of freedom of the Q statistic: 27
## P value of the Q statistic: 0.0001932118
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)   0.0337
## Intercept2: I2 (Q statistic)   0.7224
## Intercept3: I2 (Q statistic)   0.2676
## 
## Number of studies (or clusters): 10
## Number of observed statistics: 30
## Number of estimated parameters: 6
## Degrees of freedom: 24
## -2 log likelihood: -22.61046 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Extract the fixed-effects estimates
(est_fixed <- coef(random1, select="fixed"))
## Intercept1 Intercept2 Intercept3 
##  0.3777491  0.3807843  0.1704927
## Convert the estimated vector to a symmetrical matrix
## where the diagonals are fixed at 1 (for a correlation matrix)
vec2symMat(est_fixed, diag=FALSE)
##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.3777491 0.3807843
## [2,] 0.3777491 1.0000000 0.1704927
## [3,] 0.3807843 0.1704927 1.0000000

Random-effects TSSEM: Stage 2 Analysis

random2 <- tssem2(random1, Amatrix=RAM$A, Smatrix=RAM$S, Fmatrix=RAM$F)
summary(random2)
## 
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, Amatrix = Amatrix, 
##     Smatrix = Smatrix, Fmatrix = Fmatrix, diag.constraints = diag.constraints, 
##     cor.analysis = cor.analysis, intervals.type = intervals.type, 
##     mx.algebras = mx.algebras, model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##                       Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## SAT_MathONSAT_Verbal  0.325853  0.080000 0.169056 0.482649  4.0732 4.638e-05 ***
## SAT_MathONSpatial     0.322194  0.042415 0.239063 0.405325  7.5963 3.042e-14 ***
## SpatialWITHSAT_Verbal 0.170493  0.051355 0.069840 0.271146  3.3199 0.0009004 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                             Value
## Sample size                                538.00
## Chi-square of target model                   0.00
## DF of target model                           0.00
## p value of target model                      0.00
## Number of constraints imposed on "Smatrix"   0.00
## DF manually adjusted                         0.00
## Chi-square of independence model           110.81
## DF of independence model                     3.00
## RMSEA                                        0.00
## RMSEA lower 95% CI                           0.00
## RMSEA upper 95% CI                           0.00
## SRMR                                         0.00
## TLI                                          -Inf
## CFI                                          1.00
## AIC                                          0.00
## BIC                                          0.00
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## S matrix
mxEval(Smatrix, random2$mx.fit)
##           [,1]      [,2]      [,3]
## [1,] 0.7542121 0.0000000 0.0000000
## [2,] 0.0000000 1.0000000 0.1704927
## [3,] 0.0000000 0.1704927 1.0000000
## R2
mxEval(1-Smatrix, random2$mx.fit)[1,1]
## [1] 0.2457879
## Plot the model
plot(random2, color="green")

A path model for the cognitive ability to supervisor rating (1)

## Sheets including correlation matrices
sheet <- as.character(1:14)

my.R2 <- list()
for (i in sheet) {
  my.R2[[i]] <- readxl::read_excel("Hunter83.xlsx", sheet=i)
}

## Read study names and sample sizes in sheet "0"
my.study <- readxl::read_excel("Hunter83.xlsx", sheet="0")

my.R2 <- lapply(my.R2, function(x) {x <- unlist(x)
                                    x <- matrix(x, ncol=4)
                                    x <- vechs(x)
                                    vec2symMat(x, diag=FALSE)})

var.names <- c("Ability", "Knowledge", "Work sample", "Supervisor")
my.R2 <- lapply(my.R2, function(x) { dimnames(x) <- list(var.names, var.names); x})
names(my.R2) <- my.study$Study
my.R2[1:2]
## $`Campbell et al. (1973)`
##             Ability Knowledge Work sample Supervisor
## Ability        1.00      0.65        0.48       0.33
## Knowledge      0.65      1.00        0.52       0.40
## Work sample    0.48      0.52        1.00       0.23
## Supervisor     0.33      0.40        0.23       1.00
## 
## $`Corts et al. (1977)`
##             Ability Knowledge Work sample Supervisor
## Ability        1.00      0.53        0.50       0.03
## Knowledge      0.53      1.00        0.49       0.04
## Work sample    0.50      0.49        1.00       0.18
## Supervisor     0.03      0.04        0.18       1.00

A path model for cognitive ability to supervisor rating (2)

model <- "## Regression paths
          Job_knowledge ~ A2J*Ability
          Work_sample ~ A2W*Ability + J2W*Job_knowledge
          Supervisor ~ J2S*Job_knowledge + W2S*Work_sample

          ## Fix the variance of Ability at 1
          Ability ~~ 1*Ability

          ## Label the error variances of the dependent variables
          Job_knowledge ~~ Var_e_J*Job_knowledge
          Work_sample ~~ Var_e_W*Work_sample
          Supervisor ~~ Var_e_S*Supervisor"

plot(model, color="yellow", layout="spring")

RAM <- lavaan2RAM(model, obs.variables=c("Ability","Job_knowledge",
                  "Work_sample","Supervisor"))
RAM
## $A
##               Ability Job_knowledge Work_sample Supervisor
## Ability       "0"     "0"           "0"         "0"       
## Job_knowledge "0*A2J" "0"           "0"         "0"       
## Work_sample   "0*A2W" "0*J2W"       "0"         "0"       
## Supervisor    "0"     "0*J2S"       "0*W2S"     "0"       
## 
## $S
##               Ability Job_knowledge Work_sample Supervisor 
## Ability       "1"     "0"           "0"         "0"        
## Job_knowledge "0"     "0*Var_e_J"   "0"         "0"        
## Work_sample   "0"     "0"           "0*Var_e_W" "0"        
## Supervisor    "0"     "0"           "0"         "0*Var_e_S"
## 
## $F
##               Ability Job_knowledge Work_sample Supervisor
## Ability             1             0           0          0
## Job_knowledge       0             1           0          0
## Work_sample         0             0           1          0
## Supervisor          0             0           0          1
## 
## $M
##   Ability Job_knowledge Work_sample Supervisor
## 1       0             0           0          0

Random-effects TSSEM: Stage 1 Analysis

## method="REM": Random-effects model
random1 <- tssem1(Hunter83$data, Hunter83$n, method="REM", RE.type="Diag")
summary(random1)
## 
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues, 
##     "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound, 
##     I2 = I2, model.name = model.name, suppressWarnings = TRUE, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate   Std.Error      lbound      ubound z value  Pr(>|z|)    
## Intercept1  5.0045e-01  2.8249e-02  4.4508e-01  5.5581e-01 17.7157 < 2.2e-16 ***
## Intercept2  4.3124e-01  2.1560e-02  3.8898e-01  4.7350e-01 20.0020 < 2.2e-16 ***
## Intercept3  2.0615e-01  2.5700e-02  1.5578e-01  2.5652e-01  8.0214 1.110e-15 ***
## Intercept4  5.2074e-01  3.1799e-02  4.5842e-01  5.8307e-01 16.3758 < 2.2e-16 ***
## Intercept5  2.5781e-01  3.5914e-02  1.8742e-01  3.2820e-01  7.1785 7.045e-13 ***
## Intercept6  2.3471e-01  1.7505e-02  2.0040e-01  2.6902e-01 13.4085 < 2.2e-16 ***
## Tau2_1_1    6.2569e-03  3.6590e-03 -9.1449e-04  1.3428e-02  1.7100   0.08726 .  
## Tau2_2_2    2.2234e-03  2.1500e-03 -1.9905e-03  6.4373e-03  1.0342   0.30106    
## Tau2_3_3    4.3637e-03  3.0130e-03 -1.5417e-03  1.0269e-02  1.4483   0.14754    
## Tau2_4_4    7.5512e-03  4.6477e-03 -1.5581e-03  1.6661e-02  1.6247   0.10422    
## Tau2_5_5    1.0692e-02  6.0658e-03 -1.1966e-03  2.2581e-02  1.7627   0.07795 .  
## Tau2_6_6    1.0005e-10  1.7130e-03 -3.3575e-03  3.3575e-03  0.0000   1.00000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 205.225
## Degrees of freedom of the Q statistic: 60
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)   0.7209
## Intercept2: I2 (Q statistic)   0.4475
## Intercept3: I2 (Q statistic)   0.5671
## Intercept4: I2 (Q statistic)   0.7459
## Intercept5: I2 (Q statistic)   0.7691
## Intercept6: I2 (Q statistic)   0.0000
## 
## Number of studies (or clusters): 14
## Number of observed statistics: 66
## Number of estimated parameters: 12
## Degrees of freedom: 54
## -2 log likelihood: -130.2043 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Occasionally, we may encounter some error messages.
## We may try to remove the error messages by rerunning the model.
random1 <- rerun(random1)
## 
## Beginning initial fit attempt
## 
## Solution found
## 
## Solution found!  Final fit=-130.20428 (started at -130.20428)  (1 attempt(s): 1 valid, 0 errors)
summary(random1)
## 
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues, 
##     "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound, 
##     I2 = I2, model.name = model.name, suppressWarnings = TRUE, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate   Std.Error      lbound      ubound z value  Pr(>|z|)    
## Intercept1  5.0045e-01  2.8249e-02  4.4508e-01  5.5581e-01 17.7157 < 2.2e-16 ***
## Intercept2  4.3124e-01  2.1560e-02  3.8898e-01  4.7350e-01 20.0020 < 2.2e-16 ***
## Intercept3  2.0615e-01  2.5700e-02  1.5578e-01  2.5652e-01  8.0214 1.110e-15 ***
## Intercept4  5.2074e-01  3.1799e-02  4.5842e-01  5.8307e-01 16.3758 < 2.2e-16 ***
## Intercept5  2.5781e-01  3.5914e-02  1.8742e-01  3.2820e-01  7.1785 7.045e-13 ***
## Intercept6  2.3471e-01  1.7505e-02  2.0040e-01  2.6902e-01 13.4085 < 2.2e-16 ***
## Tau2_1_1    6.2569e-03  3.6590e-03 -9.1449e-04  1.3428e-02  1.7100   0.08726 .  
## Tau2_2_2    2.2234e-03  2.1500e-03 -1.9905e-03  6.4373e-03  1.0342   0.30106    
## Tau2_3_3    4.3637e-03  3.0130e-03 -1.5417e-03  1.0269e-02  1.4483   0.14754    
## Tau2_4_4    7.5512e-03  4.6477e-03 -1.5581e-03  1.6661e-02  1.6247   0.10422    
## Tau2_5_5    1.0692e-02  6.0658e-03 -1.1966e-03  2.2581e-02  1.7627   0.07795 .  
## Tau2_6_6    1.0005e-10  1.7130e-03 -3.3575e-03  3.3575e-03  0.0000   1.00000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 205.225
## Degrees of freedom of the Q statistic: 60
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)   0.7209
## Intercept2: I2 (Q statistic)   0.4475
## Intercept3: I2 (Q statistic)   0.5671
## Intercept4: I2 (Q statistic)   0.7459
## Intercept5: I2 (Q statistic)   0.7691
## Intercept6: I2 (Q statistic)   0.0000
## 
## Number of studies (or clusters): 14
## Number of observed statistics: 66
## Number of estimated parameters: 12
## Degrees of freedom: 54
## -2 log likelihood: -130.2043 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Extract the fixed-effects estimates
(est_fixed <- coef(random1, select="fixed"))
## Intercept1 Intercept2 Intercept3 Intercept4 Intercept5 Intercept6 
##  0.5004463  0.4312399  0.2061453  0.5207418  0.2578125  0.2347090
## Convert the estimated vector to a symmetrical matrix
## where the diagonals are fixed at 1 (for a correlation matrix)
vec2symMat(est_fixed, diag=FALSE)
##           [,1]      [,2]      [,3]      [,4]
## [1,] 1.0000000 0.5004463 0.4312399 0.2061453
## [2,] 0.5004463 1.0000000 0.5207418 0.2578125
## [3,] 0.4312399 0.5207418 1.0000000 0.2347090
## [4,] 0.2061453 0.2578125 0.2347090 1.0000000

Random-effects TSSEM: Stage 2 Analysis (1)

random2 <- tssem2(random1, Amatrix=RAM$A, Smatrix=RAM$S, Fmatrix=RAM$F,
                  intervals.type="LB", 
                  mx.algebras=
                  list( ind=mxAlgebra(A2J*J2S+A2J*J2W*W2S+A2W*W2S, name="ind") )  )
summary(random2)
## 
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, Amatrix = Amatrix, 
##     Smatrix = Smatrix, Fmatrix = Fmatrix, diag.constraints = diag.constraints, 
##     cor.analysis = cor.analysis, intervals.type = intervals.type, 
##     mx.algebras = mx.algebras, model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: Likelihood-based statistic
## Coefficients:
##     Estimate Std.Error   lbound   ubound z value Pr(>|z|)
## A2J 0.510783        NA 0.456554 0.564987      NA       NA
## J2S 0.224033        NA 0.136802 0.311682      NA       NA
## W2S 0.121026        NA 0.056817 0.181657      NA       NA
## A2W 0.230913        NA 0.157756 0.299775      NA       NA
## J2W 0.397800        NA 0.311509 0.484687      NA       NA
## 
## mxAlgebras objects (and their 95% likelihood-based CIs):
##             lbound  Estimate    ubound
## ind[1,1] 0.1374292 0.1669702 0.1983626
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                3975.0000
## Chi-square of target model                    3.5748
## DF of target model                            1.0000
## p value of target model                       0.0587
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model            975.2843
## DF of independence model                      6.0000
## RMSEA                                         0.0255
## RMSEA lower 95% CI                            0.0000
## RMSEA upper 95% CI                            0.0561
## SRMR                                          0.0204
## TLI                                           0.9841
## CFI                                           0.9973
## AIC                                           1.5748
## BIC                                          -4.7129
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## Display the A matrix
mxEval(Amatrix, random2$mx.fit)
##                 Ability Job_knowledge Work_sample Supervisor
## Ability       0.0000000     0.0000000   0.0000000          0
## Job_knowledge 0.5107833     0.0000000   0.0000000          0
## Work_sample   0.2309135     0.3978003   0.0000000          0
## Supervisor    0.0000000     0.2240334   0.1210258          0
## Display the S matrix
mxEval(Smatrix, random2$mx.fit)
##      [,1]      [,2]      [,3]     [,4]
## [1,]    1 0.0000000 0.0000000 0.000000
## [2,]    0 0.7391004 0.0000000 0.000000
## [3,]    0 0.0000000 0.6945954 0.000000
## [4,]    0 0.0000000 0.0000000 0.907194

Random-effects TSSEM: Stage 2 Analysis (2)

plot(random2, color="yellow", layout="spring")

Summary


  1. National Research Council (1992). Combining information: Statistical issues and opportunities for research. Washington, D.C.: National Academy Press.

  2. Brown, S. P., & Stayman, D. M. (1992). Antecedents and consequences of attitude toward the ad: A meta-analysis. Journal of Consumer Research, 19, 34-51.

  3. Premack, S. L., & Hunter, J. E. (1988). Individual unionization decisions. Psychological Bulletin, 103, 223-234.

  4. Norton, S., Cosco, T., Doyle, F., Done, J., & Sacker, A. (2013). The Hospital Anxiety and Depression Scale: A meta confirmatory factor analysis. Journal of Psychosomatic Research, 74(1), 74-81. http://doi.org/10.1016/j.jpsychores.2012.10.010

  5. Murayama, K., & Elliot, A. J. (2012). The competition-performance relation: A meta-analytic review and test of the opposing processes model of competition and performance. Psychological Bulletin, 138(6), 1035-1070. http://doi.org/10.1037/a0028324

  6. Viswesvaran, C., & Ones, D. S. (1995). Theory testing: Combining psychometric meta-analysis and structural equations modeling. Personnel Psychology, 48(4), 865-885.

  7. Becker, B. J. (1992). Using results from replicated studies to estimate linear models. Journal of Educational Statistics, 17(4), 341-362.

  8. Cheung, M. W.-L. (2014). Fixed- and random-effects meta-analytic structural equation modeling: Examples and analyses in R. Behavior Research Methods, 46(1), 29-40. http://doi.org/10.3758/s13428-013-0361-y

  9. Cheung, M. W.-L., & Chan, W. (2005). Meta-analytic structural equation modeling: A two-stage approach. Psychological Methods, 10(1), 40-64. http://doi.org/10.1037/1082-989X.10.1.40

  10. Cheung, M. W.-L. (2018). Some reflections on combining meta-analysis and structural equation modeling. Research Synthesis Methods, 0(0). https://doi.org/10.1002/jrsm.1321

  11. Cudeck, R. (1989). Analysis of correlation matrices using covariance structure models. Psychological Bulletin, 105(2), 317-327.

  12. Joreskog, K. G., & Sorbom, D. (1996). LISREL 8: A user-s reference guide. Chicago, IL: Scientific Software International, Inc.

  13. Muthen, B., Kaplan, D., & Hollis, M. (1987). On structural equation modeling with data that are not missing completely at random. Psychometrika, 52(3), 431-462. http://doi.org/10.1007/BF02294365

  14. Cheung, M. W.-L., & Chan, W. (2009). A two-stage approach to synthesizing covariance matrices in meta-analytic structural equation modeling. Structural Equation Modeling: A Multidisciplinary Journal, 16(1), 28-53. http://doi.org/10.1080/10705510802561295

  15. Bentler, P. M., & Savalei, V. (2010). Analysis of correlation structures: Current status and open problems. In S. Kolenikov, D. Steinley, & L. Thombs (Eds.), Statistics in the Social Sciences (pp. 1-36). New Jersey: John Wiley & Sons, Inc.

  16. Cheung, M. W.-L., & Chan, W. (2005). Classifying correlation matrices into relatively homogeneous subgroups: A cluster analytic approach. Educational and Psychological Measurement, 65(6), 954–979. https://doi.org/10.1177/0013164404273946

  17. Jak, S., & Cheung, M. W.-L. (2018). Testing moderator hypotheses in meta-analytic structural equation modeling using subgroup analysis. Behavior Research Methods, 50(4), 1359–1373. https://doi.org/10.3758/s13428-018-1046-3

  18. Rosenthal, R. (1991). Meta-analytic procedures for social research (Revised). Newbury Park, [Calif.]: Sage.

  19. Schmidt, F. L., & Hunter, J. E. (2015). Methods of meta-analysis: Correcting error and bias in research findings (3rd ed.). Thousand Oaks, CA: Sage.

  20. Michel, J. S., Viswesvaran, C., & Thomas, J. (2011). Conclusions from meta-analytic structural equation models generally do not change due to corrections for study artifacts. Research Synthesis Methods, 2(3), 174-187. http://doi.org/10.1002/jrsm.47

  21. McArdle, J. J., & McDonald, R. P. (1984). Some algebraic properties of the Reticular Action Model for moment structures. British Journal of Mathematical and Statistical Psychology, 37(2), 234-251. http://doi.org/10.1111/j.2044-8317.1984.tb00802.x

  22. Browne, M. W., & Cudeck, R. (1993). Alternative ways of assessing model fit. In K. A. Bollen & J. S. Long (Eds.), Testing Structural Equation Models (pp. 136-162). Newbury Park, CA: Sage.

  23. Digman, J. M. (1997). Higher-order factors of the Big Five. Journal of Personality and Social Psychology, 73(6), 1246-1256. http://doi.org/10.1037/0022-3514.73.6.1246

  24. Epskamp, S (2017). semPlot: Path diagrams and visual analysis of various SEM packages’ output. R package version 1.1.0. http://CRAN.R-project.org/package=semPlot

  25. Becker, B. J., & Schram, C. M. (1994). Examining explanatory models through research synthesis. In H. Cooper & L. V. Hedges (Eds.), The handbook of research synthesis (pp. 357-381). New York: Russell Sage Foundation.